Script projet CIV6708 H23

Author

Victor Nunzi

Published

March 18, 2023

Script Quarto ayant pour but de documenter les traitements effectués dans le cadre du projet CIV6708.

Date du dernier render : 2023-03-24

1 Librairies

library(sf)
library(tidyverse)
library(stringi)
library(DT)
library(ggplot2)
library(tmap)

2 Recensement

2.1 Données

2.1.1 Population et tranches d’âge

Données brut du recensement avec la population par tranches d’âge, par Aires de diffusion (AD).

Source : recensement 2021

recensement.age <- read.csv("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/StatCan/recensement_age_AD.csv", 
                            sep = ";", fileEncoding = "UTF-8-BOM")

Echelle : Canada entier (907752 lignes).

Affichage des premières lignes :

2.1.2 Limite des AD

On charge les limites géographiques des AD de Drummondville :

sf.drummond.ville.AD <- st_read("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/SDR_drummondville/AD_drummondville.shp")
Reading layer `AD_drummondville' from data source 
  `C:\Users\victo\OneDrive - polymtl.ca\CIV6708\PROJET\SDR_drummondville\AD_drummondville.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 137 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 7684020 ymin: 1299056 xmax: 7719001 ymax: 1324995
Projected CRS: NAD83 / Statistics Canada Lambert

Affichage des éléments :

2.2 Traitement

On commence par formater les noms des colonnes :

colnames <- colnames(recensement.age)

colnames <- stri_trans_general(colnames, id = "Latin-ASCII") # retire les accents 
colnames <- tolower(colnames)                                # met tout en minuscules
colnames <- str_replace_all(colnames, "\\.", "_")            # remplace les points par des underscores

colnames(recensement.age) <- colnames

On transforme les données :

recensement.age <- recensement.age %>%
  
  select( # sélection des colonnes intéressantes
    c(
      geo,
      dguid,
      age__16_,
      valeur
    )
  ) %>%
  
  filter( # sélection des données de Drummondville seulement
    geo %in% sf.drummond.ville.AD$ADIDU
  ) %>%
  
  pivot_wider( # on pivote les données dans la largeur (afin d'avoir une ligne par AD)
    names_from = "age__16_", values_from = "valeur"
    )

On procède aux mêmes traitements que précédemment sur les nouvelles colonnes :

colnames <- colnames(recensement.age)

colnames <- stri_trans_general(colnames, id = "Latin-ASCII")
colnames <- tolower(colnames)
colnames <- str_replace_all(colnames, "\\.", "_")
colnames <- str_replace_all(colnames, " ", "_")
colnames <- str_remove(colnames, "-")

colnames(recensement.age) <- colnames

On remplace les NA par des 0 :

recensement.age[is.na(recensement.age)] <- 0

On affiche les données traitées :

2.3 Enregistrement

write.csv(recensement.age, "C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/StatCan/recensement_age_AD_tri.csv", row.names = F)

3 Association de la population aux bâtiments

3.1 Données

On charge les données de population par tranches d’âge (traitées précédemment) :

recensement.age <- read.csv("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/StatCan/recensement_age_AD_tri.csv")

On charge également les bâtiments résidentiels, extraits de OSM.

Clé utilisée : building

Valeurs : apartments, cabin, detached, dormitory, semidetached_house, static_caravan

sf.accomodations <- st_read("C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/OSM/buildings_accomodations_drummondville.shp")
Reading layer `buildings_accomodations_drummondville' from data source 
  `C:\Users\victo\OneDrive - polymtl.ca\CIV6708\PROJET\OSM\buildings_accomodations_drummondville.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 5329 features and 31 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -72.64695 ymin: 45.75518 xmax: -72.28676 ymax: 46.02368
Geodetic CRS:  WGS 84

Avant d’afficher les données, on fait quelques sélections et changements dans les noms des champs :

sf.accomodations <- sf.accomodations %>%
  select( # sélection des champs qui nous intéressent
    c(
      full_id:building,
      building_1:geometry
    )
  ) %>%
  
  rename(
    flats = building_1 # on corrige l'abbréviation automatique du format shapefile
  )

Affichage :

Pour information, le champ acc_number a été créé dans QGIS, et vaut 1 si flats est vide, 0 sinon.

3.2 Traitement

On commence par changer le CRS des données OSM, pour concorder avec celui des données de StatCan :

sf.accomodations <- st_transform(sf.accomodations, crs = 3347)

On associe à chaque bâtiment le code de l’AD dans laquelle son centroïde se trouve :

sf.accomodations.points <- st_join(st_centroid(sf.accomodations), sf.drummond.ville.AD[c("IDUGD")],
                                   join = st_within)
Warning: st_centroid assumes attributes are constant over geometries
sf.accomodations$IDUGD <- 
  sf.accomodations.points$IDUGD[match(sf.accomodations$full_id,
                                      sf.accomodations.points$full_id)]

On ajoute à chaque bâtiment les données de population pour l’AD dans laquelle il se trouve :

sf.accomodations <- left_join(sf.accomodations, recensement.age,
                              by = c("IDUGD" = "dguid"))

Affichage des premières lignes de la table des bâtiments :

On a actuellement, pour chaque bâtiment, les données de population pour l’AD au complet. On souhaite répartir ces données à l’échelle de chaque bâtiment, en se basant sur deux conditions :

  • les sommes des données de population pour tous les bâtiments d’une AD doit être égale aux données de StatCan ;
  • le nombre de personnes attribuées à chaque bâtiment doit être proportionnel au nombre de logements de ce bâtiments.

On va pour cela utiliser la fonction suivante, qui permet de générer un vecteur de nombre entier à partir de la somme de ce vecteur et d’un vecteur de probabilités.

Source : https://stackoverflow.com/questions/57278562/generating-randomly-a-vector-of-integers-that-sum-up-to-n-with-a-given-vector-of

rand_vect <- function(N, M, sd = 1, pos.only = TRUE,prob) {
  vec <- rbinom(N,1,prob)
  deviation <- M - sum(vec)
  for (. in seq_len(abs(deviation))) {
    vec[i] <- vec[i <- sample(N, 1,replace=TRUE,prob)] + sign(deviation)
  }
  if (pos.only) while (any(vec < 0)) {
    negs <- vec < 0
    pos  <- vec > 0
    vec[negs][i] <- vec[negs][i <- sample(sum(negs), 1)] + 1
    vec[pos][i]  <- vec[pos ][i <- sample(sum(pos ), 1)] - 1
  }
  vec
}

On attribue les personnes à chaque bâtiment :

df.accomodations <- st_drop_geometry(sf.accomodations) %>% # on travaille avec un dataframe plutôt qu'un sf
  
  group_by(IDUGD) %>% # on travaille par AD
  
  mutate(
    nb_bat_AD    = n(),                                  # nombre de bâtiments dans l'AD
    share_log_AD = round(acc_number/sum(acc_number), 3), # part des logements de l'AD représentée par le bâtiment
    total_pop_AD = total__age,
    
    across(X15_a_19_ans:X85_ans_et_plus,    # opération effectuée à l'identique pour chaque tranche d'âge
           ~ rand_vect(N = n(),             # nombre de bâtiments à populer dans l'AD
                       M = unique(.x),      # nombre de personne dans l'AD pour cette tranche d'âge
                       prob = share_log_AD) # proportion des logements représentée par chaque bâtiment pour l'AD
           )
  ) %>%
  
  rowwise() %>% mutate(
    total__age = sum(c_across(X15_a_19_ans:X85_ans_et_plus)) # pour le total, on somme ce qu'on a obtenu pour chaque tranche d'âges
  ) %>%
  
  # quelques opérations pour réorganiser les colonnes :
  select(-geo) %>%
  
  relocate(IDUGD) %>%
  
  relocate(nb_bat_AD, 
           share_log_AD,
           total_pop_AD,
           .before = total__age)

sf.accomodations <- left_join(sf.accomodations[c("full_id")], df.accomodations) %>%
  relocate(IDUGD) # on ajoute les traitements aux géométries d'origine
Joining with `by = join_by(full_id)`
rm(df.accomodations)

Affichage de la table complète :

Warning in instance$preRenderHook(instance): It seems your data is too big for
client-side DataTables. You may consider server-side processing:
https://rstudio.github.io/DT/server.html

Carte :

La suite des traitements (production d’une carte de chaleur) va être effectuée dans QGIS.

3.3 Enregistrement

st_write(sf.accomodations, "C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/OSM/accomodations_drummondville.shp",
         delete_layer = T)
Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
Shapefile driver
Deleting layer `accomodations_drummondville' using driver `ESRI Shapefile'
Writing layer `accomodations_drummondville' to data source 
  `C:/Users/victo/OneDrive - polymtl.ca/CIV6708/PROJET/OSM/accomodations_drummondville.shp' using driver `ESRI Shapefile'
Writing 5329 features with 26 fields and geometry type Multi Polygon.